<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of demse2</title>
  <meta name="keywords" content="demse2">
  <meta name="description" content="DEMSE2  Demonstrate state estimation on a simple scalar nonlinear (time variant) problem">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">examples</a> &gt; <a href="#">state_estimation</a> &gt; demse2.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../../menu.html"><img alt="<" border="0" src="../../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\examples\state_estimation&nbsp;<img alt=">" border="0" src="../../../right.png"></a></td></tr></table>-->

<h1>demse2
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>DEMSE2  Demonstrate state estimation on a simple scalar nonlinear (time variant) problem</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> DEMSE2  Demonstrate state estimation on a simple scalar nonlinear (time variant) problem

   See also
   GSSM_N1
   Copyright (c) Oregon Health &amp; Science University (2006)

   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for
   academic use only (see included license file) and can be obtained from
   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the
   software should contact rebel@csee.ogi.edu for commercial licensing information.

   See LICENSE (which should be part of the main toolkit distribution) for more
   detail.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>	ADDRELPATH  Add a relative path which gets expanded into a absolute path</li><li><a href="../../.././ReBEL-0.2.7/core/cdkf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = cdkf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">cdkf</a>	CDKF  Central Difference Kalman Filter  (Sigma-Point Kalman Filter variant)</li><li><a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>	CVECREP  Column vector replicate</li><li><a href="../../.././ReBEL-0.2.7/core/ekf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ekf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ekf</a>	EKF  Extended Kalman Filter</li><li><a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>	GENINFDS  Generate inference data structure from a generalized state space model and user defined inference parameters.</li><li><a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>	GENNOISEDS    Generates a NoiseDS data structure describing a noise source.</li><li><a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>	GENSYSNOISEDS  Generate process and observation noise data structures for a given InferenceDS data structure</li><li><a href="../../.././ReBEL-0.2.7/core/gmmfit.html" class="code" title="function [gmmDS, leb] = gmmfit(X, M, tt, cov_type, check_cov, display, W)">gmmfit</a>	GMMFIT   Fit a Gaussian mixture model (GMM) with M components to dataset X</li><li><a href="../../.././ReBEL-0.2.7/core/gmsppf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise, extra] = gmsppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">gmsppf</a>	GMSPPF  Gaussian Mixture Sigma-Point Particle Filter</li><li><a href="../../.././ReBEL-0.2.7/core/gspf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise] = gspf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">gspf</a>	GSPF  Gaussian Sum Particle Filter</li><li><a href="../../.././ReBEL-0.2.7/core/pf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise] = pf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">pf</a>	PF  Generic Particle Filter</li><li><a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>	REMRELPATH  Remove a relative path (which gets expanded into a absolute path)</li><li><a href="../../.././ReBEL-0.2.7/core/sppf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise] = sppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">sppf</a>	SPPF  Sigma-Point Particle Filter.</li><li><a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>	SRCDKF  Square Root Central Difference Kalman Filter (Sigma-Point Kalman Filter variant)</li><li><a href="../../.././ReBEL-0.2.7/core/srukf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srukf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srukf</a>	SRUKF  Square Root Unscented Kalman Filter (Sigma-Point Kalman Filter variant)</li><li><a href="../../.././ReBEL-0.2.7/core/ukf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ukf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ukf</a>	UKF  Unscented Kalman Filter</li><li><a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_n1.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_n1</a>	GSSM_N1  Generalized state space model for simple nonlinear system</li></ul>
This function is called by:
<ul style="list-style-image:url(../../../matlabicon.gif)">
</ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">% DEMSE2  Demonstrate state estimation on a simple scalar nonlinear (time variant) problem</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%   See also</span>
0004 <span class="comment">%   GSSM_N1</span>
0005 <span class="comment">%   Copyright (c) Oregon Health &amp; Science University (2006)</span>
0006 <span class="comment">%</span>
0007 <span class="comment">%   This file is part of the ReBEL Toolkit. The ReBEL Toolkit is available free for</span>
0008 <span class="comment">%   academic use only (see included license file) and can be obtained from</span>
0009 <span class="comment">%   http://choosh.csee.ogi.edu/rebel/.  Businesses wishing to obtain a copy of the</span>
0010 <span class="comment">%   software should contact rebel@csee.ogi.edu for commercial licensing information.</span>
0011 <span class="comment">%</span>
0012 <span class="comment">%   See LICENSE (which should be part of the main toolkit distribution) for more</span>
0013 <span class="comment">%   detail.</span>
0014 
0015 <span class="comment">%=============================================================================================</span>
0016 
0017 clc;
0018 clear all;
0019 
0020 fprintf(<span class="string">'\nDEMSE2 : This demonstration shows how the ReBEL toolkit is used for simple state estimation\n'</span>);
0021 fprintf(<span class="string">'         on a scalar nonlinear problem.\n\n'</span>);
0022 
0023 
0024 <span class="comment">%--- General setup</span>
0025 
0026 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../gssm'</span>);         <span class="comment">% add relative search path to example GSSM files to MATLABPATH</span>
0027 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../data'</span>);         <span class="comment">% add relative search path to example data files to MATLABPATH</span>
0028 <a href="../../.././ReBEL-0.2.7/core/addrelpath.html" class="code" title="function addrelpath(path_string)">addrelpath</a>(<span class="string">'../../netlab'</span>);    <span class="comment">% Some of the algorithms requires Netlab functions</span>
0029 
0030 <span class="comment">%--- Ask the user which inference algorithm to use</span>
0031 <span class="comment">%ftype = input('Type of estimator [ ekf, ukf, cdkf, srukf, srcdkf, pf, gspf, gmsppf or sppf ] ? ','s');</span>
0032 <span class="comment">%if ~stringmatch(ftype,{'ekf','cdkf','ukf','srukf','srcdkf','pf','gspf','sppf','gmsppf'})</span>
0033 <span class="comment">%    error('That estimator/filter type is not recognized.');</span>
0034 <span class="comment">%end</span>
0035 
0036 <span class="comment">%--- Compare these algorithms...</span>
0037 lftype={<span class="string">'ekf'</span>,<span class="string">'cdkf'</span>,<span class="string">'ukf'</span>,<span class="string">'srukf'</span>,<span class="string">'srcdkf'</span>,<span class="string">'pf'</span>,<span class="string">'gspf'</span>,<span class="string">'sppf'</span>,<span class="string">'gmsppf'</span>};
0038 
0039 number_of_runs = input(<span class="string">'Number of independent runs ? '</span>);
0040 
0041 mean_RMSE = zeros(1,9);   <span class="comment">% buffer for MC results for each algorithm</span>
0042 var_RMSE  = zeros(1,9);   <span class="comment">%     &quot;                              &quot;</span>
0043 
0044 <span class="keyword">for</span> jj=1:9,
0045 
0046   ftype=lftype{jj};
0047 
0048   disp([<span class="string">'['</span> ftype <span class="string">']'</span>]);
0049 
0050   <span class="comment">%--- Initialise GSSM model from external system description script.</span>
0051 
0052   model = <a href="../../.././ReBEL-0.2.7/examples/gssm/gssm_n1.html" class="code" title="function [varargout] = model_interface(func, varargin)">gssm_n1</a>(<span class="string">'init'</span>);
0053 
0054   Arg.type = <span class="string">'state'</span>;                                  <span class="comment">% inference type (state estimation)</span>
0055   Arg.tag = <span class="string">'State estimation for GSSM_N1 system.'</span>;    <span class="comment">% arbitrary ID tag</span>
0056   Arg.model = model;                                   <span class="comment">% GSSM data structure of external system</span>
0057   Arg.algorithm = ftype;                               <span class="comment">% set inference algorithm to be used</span>
0058 
0059   InfDS = <a href="../../.././ReBEL-0.2.7/core/geninfds.html" class="code" title="function InferenceDS = geninfds(ArgDS)">geninfds</a>(Arg);                               <span class="comment">% Create inference data structure and</span>
0060   [pNoise, oNoise, InfDS] = <a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>(InfDS,ftype);       <span class="comment">% generate process and observation noise sources</span>
0061 
0062 
0063   <span class="comment">%--- Loop over number of independent runs</span>
0064 
0065   <span class="keyword">for</span> k=1:number_of_runs
0066 
0067     randn(<span class="string">'state'</span>,sum(100*clock));          <span class="comment">% stir the pot... shuffle the deck :-)</span>
0068     rand(<span class="string">'state'</span>,sum(100*clock));
0069 
0070     <span class="comment">%--- Generate some data</span>
0071 
0072     N  = 60;                                                <span class="comment">% number of datapoints</span>
0073     X  = zeros(model.statedim,N);                           <span class="comment">% state data buffer</span>
0074     y  = zeros(model.obsdim,N);                             <span class="comment">% observation data buffer</span>
0075 
0076     pnoise = model.pNoise.sample( model.pNoise, N);   <span class="comment">% generate process noise</span>
0077     onoise = model.oNoise.sample( model.oNoise, N);   <span class="comment">% generate observation noise</span>
0078 
0079     X(1) = 1;                                               <span class="comment">% initial state</span>
0080     y(1) = model.hfun( model, X(1), onoise(1), 1);    <span class="comment">% observation of initial state</span>
0081     <span class="keyword">for</span> j=2:N,
0082       X(j) = model.ffun( model, X(:,j-1), pnoise(j-1), j-1);
0083       y(j) = model.hfun( model, X(:,j), onoise(j), j);
0084     <span class="keyword">end</span>
0085 
0086     U1 = [0:N-1];
0087     U2 = [1:N];
0088 
0089     <span class="comment">%--- Setup runtime buffers</span>
0090 
0091     Xh = zeros(InfDS.statedim,N);          <span class="comment">% state estimation buffer</span>
0092     Xh(:,1) = 1;                           <span class="comment">% initial estimate of state E[X(0)]</span>
0093     Px = 3/4*eye(InfDS.statedim);          <span class="comment">% initial state covariance</span>
0094 
0095     <span class="comment">%--- Call inference algorithm / estimator</span>
0096 
0097     <span class="keyword">switch</span> ftype
0098 
0099 
0100       <span class="comment">%------------------- Extended Kalman Filter ------------------------------------</span>
0101       <span class="keyword">case</span> <span class="string">'ekf'</span>
0102 
0103         [Xh, Px] = <a href="../../.././ReBEL-0.2.7/core/ekf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ekf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ekf</a>(Xh(:,1), Px, pNoise, oNoise, y, U1, U2, InfDS);
0104 
0105 
0106       <span class="comment">%------------------- Unscented Kalman Filter -----------------------------------</span>
0107       <span class="keyword">case</span> <span class="string">'ukf'</span>
0108 
0109         alpha = 1;         <span class="comment">% scale factor (UKF parameter)</span>
0110         beta  = 2;         <span class="comment">% optimal setting for Gaussian priors (UKF parameter)</span>
0111         kappa = 0;         <span class="comment">% optimal for state dimension=2 (UKF parameter)</span>
0112 
0113         InfDS.spkfParams = [alpha beta kappa];
0114 
0115         [Xh, Px] = <a href="../../.././ReBEL-0.2.7/core/ukf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = ukf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">ukf</a>(Xh(:,1), Px, pNoise, oNoise, y, U1, U2, InfDS);
0116 
0117 
0118       <span class="comment">%------------------- Central Difference Kalman Filter ---------------------------</span>
0119       <span class="keyword">case</span> <span class="string">'cdkf'</span>
0120 
0121         InfDS.spkfParams = sqrt(3);    <span class="comment">% scale factor (CDKF parameter h)</span>
0122 
0123         [Xh, Px] = <a href="../../.././ReBEL-0.2.7/core/cdkf.html" class="code" title="function [xh, Px, pNoise, oNoise, InternalVariablesDS] = cdkf(state, Pstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">cdkf</a>(Xh(:,1), Px, pNoise, oNoise, y, U1, U2, InfDS);
0124 
0125 
0126       <span class="comment">%------------------- Square Root Unscented Kalman Filter ------------------------</span>
0127       <span class="keyword">case</span> <span class="string">'srukf'</span>
0128 
0129         alpha = 1;         <span class="comment">% scale factor (UKF parameter)</span>
0130         beta  = 2;         <span class="comment">% optimal setting for Gaussian priors (UKF parameter)</span>
0131         kappa = 0;         <span class="comment">% optimal for state dimension=2 (UKF parameter)</span>
0132 
0133         Sx = chol(Px)';
0134 
0135         InfDS.spkfParams = [alpha beta kappa];
0136 
0137         [Xh, Sx] = <a href="../../.././ReBEL-0.2.7/core/srukf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srukf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srukf</a>(Xh(:,1), Sx, pNoise, oNoise, y, U1, U2, InfDS);
0138 
0139 
0140       <span class="comment">%------------------- Square Root Central Difference Kalman Filter ---------------</span>
0141       <span class="keyword">case</span> <span class="string">'srcdkf'</span>
0142 
0143         InfDS.spkfParams  = sqrt(3);    <span class="comment">% scale factor (CDKF parameter h)</span>
0144 
0145         Sx = chol(Px)';
0146 
0147         [Xh, Sx] = <a href="../../.././ReBEL-0.2.7/core/srcdkf.html" class="code" title="function [xh, Sx, pNoise, oNoise, InternalVariablesDS] = srcdkf(state, Sstate, pNoise, oNoise, obs, U1, U2, InferenceDS)">srcdkf</a>(Xh(:,1), Sx, pNoise, oNoise, y, U1, U2, InfDS);
0148 
0149 
0150       <span class="comment">%------------------- Generic Particle Filter (a.k.a Bootstrap-filter of CONDENSATION -----------</span>
0151       <span class="keyword">case</span> <span class="string">'pf'</span>
0152 
0153         M = 200;                             <span class="comment">% number of particles</span>
0154         ParticleFiltDS.N = M;
0155         ParticleFiltDS.particles = randn(InfDS.statedim,M)+<a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(Xh(:,1),M);  <span class="comment">% initialize particles</span>
0156         ParticleFiltDS.weights = <a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(1/M,M); <span class="comment">% initialize weights</span>
0157 
0158         InfDS.resampleThreshold = 0.5;    <span class="comment">% set resample threshold</span>
0159         InfDS.estimateType = <span class="string">'mean'</span>;      <span class="comment">% estimate type for Xh</span>
0160 
0161         [Xh, ParticleFiltDS] = <a href="../../.././ReBEL-0.2.7/core/pf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise] = pf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">pf</a>(ParticleFiltDS, pNoise, oNoise, y, U1, U2, InfDS);
0162 
0163       <span class="comment">%------------------- Gaussian-Sum Particle Filter ---------------------------------------------</span>
0164       <span class="keyword">case</span> <span class="string">'gspf'</span>
0165 
0166         M = 200;                             <span class="comment">% number of particles</span>
0167         ParticleFiltDS.N = M;
0168 
0169         initialParticles = randn(InfDS.statedim,M)+<a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(Xh(:,1),M);  <span class="comment">% initialize particles</span>
0170 
0171         ParticleFiltDS.stateGMM = <a href="../../.././ReBEL-0.2.7/core/gmmfit.html" class="code" title="function [gmmDS, leb] = gmmfit(X, M, tt, cov_type, check_cov, display, W)">gmmfit</a>(initialParticles, 2, [0.001 10], <span class="string">'sqrt'</span>);  <span class="comment">% fit a 3 component GMM to initial state distribution</span>
0172 
0173         InfDS.estimateType = <span class="string">'mean'</span>;      <span class="comment">% estimate type for Xh</span>
0174         InfDS.threshold = 0.001;
0175 
0176         Arg.type=<span class="string">'gmm'</span>;
0177         Arg.cov_type=<span class="string">'sqrt'</span>;
0178         Arg.dim=model.Vdim;
0179         Arg.M = 2;
0180         Arg.mu = <a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(model.pNoise.mu,Arg.M);
0181         Arg.cov = zeros(Arg.dim,Arg.dim,Arg.M);
0182         Arg.cov(:,:,1) = 2*model.pNoise.cov(:,:,1);
0183         Arg.cov(:,:,2) = 0.5*model.pNoise.cov(:,:,1);
0184         Arg.weights = [0.5 0.5];
0185         pNoise = <a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>(Arg);
0186 
0187         [Xh, ParticleFiltDS] = <a href="../../.././ReBEL-0.2.7/core/gspf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise] = gspf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">gspf</a>(ParticleFiltDS, pNoise, oNoise, y, U1, U2, InfDS);
0188 
0189       <span class="comment">%------------------- Sigma-Point Bayes Filter ---------------------------------------------</span>
0190       <span class="keyword">case</span> <span class="string">'gmsppf'</span>
0191 
0192         M = 200;
0193         ParticleFiltDS.N = M;            <span class="comment">% number of particles</span>
0194 
0195         initialParticles = randn(InfDS.statedim,M)+<a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(Xh(:,1),M);  <span class="comment">% initialize particles</span>
0196 
0197         tempCov = zeros(1,1,2); tempCov(:,:,1) = sqrt(2); tempCov(:,:,2)=1;
0198 
0199         ParticleFiltDS.stateGMM = <a href="../../.././ReBEL-0.2.7/core/gmmfit.html" class="code" title="function [gmmDS, leb] = gmmfit(X, M, tt, cov_type, check_cov, display, W)">gmmfit</a>(initialParticles, 3, [0.001 10], <span class="string">'sqrt'</span>);  <span class="comment">% fit a 3 component GMM to initial state distribution</span>
0200 
0201         InfDS.estimateType = <span class="string">'mean'</span>;    <span class="comment">% estimate type for Xh</span>
0202 
0203         InfDS.spkfType = <span class="string">'srcdkf'</span>;      <span class="comment">% Type of SPKF to use inside SPPF (note that ParticleFiltDS.particlesCov should comply)</span>
0204         InfDS.spkfParams  = sqrt(3);    <span class="comment">% scale factor (CDKF parameter h)</span>
0205 
0206         Arg.type=<span class="string">'gmm'</span>;
0207         Arg.cov_type=<span class="string">'sqrt'</span>;
0208         Arg.dim=model.Vdim;
0209         Arg.M = 2;
0210         Arg.mu = <a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(model.pNoise.mu,Arg.M);
0211         Arg.cov = zeros(Arg.dim,Arg.dim,Arg.M);
0212         Arg.cov(:,:,1) = 2*model.pNoise.cov(:,:,1);
0213         Arg.cov(:,:,2) = 0.5*model.pNoise.cov(:,:,1);
0214         Arg.weights = [0.5 0.5];
0215         pNoise = <a href="../../.././ReBEL-0.2.7/core/gennoiseds.html" class="code" title="function NoiseDS = gennoiseds(ArgDS)">gennoiseds</a>(Arg);
0216 
0217         [Xh, ParticleFiltDS] = <a href="../../.././ReBEL-0.2.7/core/gmsppf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise, extra] = gmsppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">gmsppf</a>(ParticleFiltDS, pNoise, oNoise, y, U1, U2, InfDS);
0218 
0219 
0220 
0221       <span class="comment">%------------------- Sigma-Point Particle Filter -----------------------------------------------</span>
0222       <span class="keyword">case</span> <span class="string">'sppf'</span>
0223 
0224         M = 200;                             <span class="comment">% number of particles</span>
0225         ParticleFiltDS.N = M;
0226         ParticleFiltDS.particles  = <a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(Xh(:,1),M);  <span class="comment">% initialize particle means</span>
0227         ParticleFiltDS.particlesCov = repmat(eye(InfDS.statedim),[1 1 M]);       <span class="comment">% particle covariances</span>
0228 
0229         pNoiseGAUS.cov = sqrt(2*3/4);
0230         oNoiseGAUS.cov = sqrt(1e-1);
0231 
0232         [pNoiseGAUS, oNoiseGAUS, foo] = <a href="../../.././ReBEL-0.2.7/core/gensysnoiseds.html" class="code" title="function [pNoise, oNoise, InferenceDS] = gensysnoiseds(InferenceDS, estimatorType, pNoiseAdaptMethod, pNoiseAdaptParams,oNoiseAdaptMethod, oNoiseAdaptParams)">gensysnoiseds</a>(InfDS,<span class="string">'srukf'</span>);
0233 
0234         ParticleFiltDS.pNoise = pNoiseGAUS;
0235         ParticleFiltDS.oNoise = oNoiseGAUS;
0236         ParticleFiltDS.weights = <a href="../../.././ReBEL-0.2.7/core/cvecrep.html" class="code" title="function m = cvecrep(v,c)">cvecrep</a>(1/M,M); <span class="comment">% initialize weights</span>
0237 
0238         InfDS.spkfType = <span class="string">'srukf'</span>;         <span class="comment">% Type of SPKF to use (note that ParticleFiltDS.particlesP should comply)</span>
0239         <span class="comment">%InfDS.spkfParams = [sqrt(3)];</span>
0240         InfDS.spkfParams = [1 0 2];
0241         InfDS.resampleThreshold = 1;    <span class="comment">% set resample threshold</span>
0242         InfDS.estimateType = <span class="string">'mean'</span>;      <span class="comment">% estimate type for Xh</span>
0243 
0244         [Xh, ParticleFiltDS] = <a href="../../.././ReBEL-0.2.7/core/sppf.html" class="code" title="function [estimate, ParticleFilterDS, pNoise, oNoise] = sppf(ParticleFilterDS, pNoise, oNoise, obs, U1, U2, InferenceDS)">sppf</a>(ParticleFiltDS, pNoise, oNoise, y, U1, U2, InfDS);
0245 
0246     <span class="keyword">end</span>
0247 
0248     <span class="comment">%--- Plot results</span>
0249 
0250     figure(1); clf;
0251     p1 = plot(X(1,:)); hold on
0252     p2 = plot(y,<span class="string">'g+'</span>);
0253     p3 = plot(Xh(1,:),<span class="string">'r'</span>); hold off;
0254     legend([p1 p2 p3],<span class="string">'clean'</span>,<span class="string">'noisy'</span>,[ftype <span class="string">' estimate'</span>]);
0255     xlabel(<span class="string">'time'</span>);
0256     title(<span class="string">'DEMSE2 : Nonlinear Time Variant State Estimation (non Gaussian noise)'</span>);
0257 
0258     drawnow
0259 
0260     <span class="comment">%--- Calculate mean square estimation error</span>
0261 
0262     rmse(k) = sqrt(mean((Xh(1,2:end)-X(1,2:end)).^2));
0263     fprintf(<span class="string">'%d:%d  Root-mean-square-error (RMSE) of estimate : %4.3f\n'</span>, k, number_of_runs, rmse(k));
0264 
0265   <span class="keyword">end</span>
0266 
0267   disp(<span class="string">' '</span>);
0268   mean_RMSE(jj) = mean(rmse);
0269   var_RMSE(jj)  = var(rmse);
0270 
0271 <span class="keyword">end</span>
0272 
0273 
0274 <span class="comment">%--- Summary of results</span>
0275 disp(<span class="string">' '</span>);
0276 disp(<span class="string">' '</span>);
0277 <span class="keyword">for</span> jj=1:9,
0278    disp([lftype{jj} <span class="string">' - RMSE (mean) : '</span> num2str(mean_RMSE(jj)) <span class="string">'   RMSE (var) : '</span> num2str(var_RMSE(jj))]);
0279 <span class="keyword">end</span>
0280 
0281 <span class="comment">%--- House keeping</span>
0282 
0283 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../gssm'</span>);       <span class="comment">% remove relative search path to example GSSM files from MATLABPATH</span>
0284 <a href="../../.././ReBEL-0.2.7/core/remrelpath.html" class="code" title="function remrelpath(path_string)">remrelpath</a>(<span class="string">'../data'</span>);       <span class="comment">% remove relative search path to example data files from MATLABPATH</span></pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>